---
title: Multidimensional arrays
---

Common Lisp has native support for multidimensional arrays, with some
special treatment for 1-D arrays, called `vectors`. Arrays can be
*generalised* and contain any type (`element-type t`), or they
can be *specialised* to contain specific types such as `single-float`
or `integer`. A good place to start is
[Practical Common Lisp Chapter 11, Collections](http://www.gigamonkeys.com/book/collections.html) by
Peter Seibel.

A quick reference to some common operations on arrays is given in the
section on [Arrays and vectors](data-structures.html).

Some libraries available on [Quicklisp](https://www.quicklisp.org/beta/) for manipulating arrays:

- [array-operations](https://github.com/Lisp-Stat/array-operations) maintained by @Symbolics defines
  functions `generate`, `permute`, `displace`, `flatten`, `split`,
  `combine`, `reshape`. It also defines `each`, for element-wise
  operations. This is a fork of [bendudson/array-operations](https://github.com/bendudson/array-operations) which is a fork of
  [tpapp/array-operations](https://github.com/tpapp/array-operations), the original
  author.
- [cmu-infix](https://github.com/rigetticomputing/cmu-infix) includes
  array indexing syntax for multidimensional arrays.
- [lla](https://github.com/tpapp/lla) is a library for linear algebra, calling BLAS and LAPACK
  libraries. It differs from most CL linear algebra packages in using
  intuitive function names, and can operate on native arrays as well as
  CLOS objects.

This page covers what can be done with the built-in multidimensional
arrays, but there are limitations. In particular:

* Interoperability with foreign language arrays, for example when
  calling libraries such as BLAS, LAPACK or GSL.
* Extending arithmetic and other mathematical operators to handle
  arrays, for example so that `(+ a b)` works
  when `a` and/or `b` are arrays.

Both of these problems can be solved by using CLOS to define an
extended array class, with native arrays as a special case.
Some libraries available through
[quicklisp](https://www.quicklisp.org/beta/) which take this approach
are:

* [matlisp](https://github.com/bharath1097/matlisp/), some of which is
  described in sections below.
* [MGL-MAT](https://github.com/melisgl/mgl-mat), which has a manual
  and provides bindings to BLAS and CUDA. This is used in a machine
  learning library [MGL](https://github.com/melisgl/mgl).
* [cl-ana](https://github.com/ghollisjr/cl-ana/wiki), a data analysis
  package with a manual, which includes operations on arrays.
* [Antik](https://www.common-lisp.net/project/antik/), used in
  [GSLL](https://common-lisp.net/project/gsll/), a binding to the GNU
  Scientific Library.

A relatively new but actively developed package is
[MAGICL](https://github.com/rigetticomputing/magicl), which provides
wrappers around BLAS and LAPACK libraries. At the time of writing this
package is not on Quicklisp, and only works under SBCL and CCL. It
seems to be particularly focused on complex arrays, but not
exclusively.
To install, clone the repository in your quicklisp `local-projects`
directory e.g. under Linux/Unix:

~~~bash
$ cd ~/quicklisp/local-projects
$ git clone https://github.com/rigetticomputing/magicl.git
~~~

Instructions for installing dependencies (BLAS, LAPACK and Expokit)
are given on the [github web pages](https://github.com/rigetticomputing/magicl).
Low-level routines wrap foreign functions, so have the Fortran names
e.g `magicl.lapack-cffi::%zgetrf`. Higher-level interfaces to some of
these functions also exist, see the
[source directory](https://github.com/rigetti/magicl/blob/master/src/high-level/) and [documentation](https://github.com/quil-lang/magicl/blob/master/doc/high-level.md).

Taking this further, domain specific languages have been built on Common
Lisp, which can be used for numerical calculations with arrays.
At the time of writing the most widely used and supported of these are:

* [Maxima](http://maxima.sourceforge.net/documentation.html)
* [Axiom](https://github.com/daly/axiom)


[CLASP](https://github.com/drmeister/clasp) is a project
which aims to ease interoperability of Common Lisp with other
languages (particularly C++), by using [LLVM](http://llvm.org/).
One of the main applications of this project is to numerical/scientific
computing.

## Creating

The function [CLHS: make-array](http://clhs.lisp.se/Body/f_mk_ar.htm)
can create arrays filled with a single value

~~~lisp
* (defparameter *my-array* (make-array '(3 2) :initial-element 1.0))
*MY-ARRAY*
* *my-array*
#2A((1.0 1.0) (1.0 1.0) (1.0 1.0))
~~~

More complicated array values can be generated by first making an
array, and then iterating over the elements to fill in the values (see
section below on element access).

The [array-operations](https://github.com/tpapp/array-operations)
library provides `generate`, a convenient
function for creating arrays which wraps this iteration.

~~~lisp
* (ql:quickload :array-operations)
To load "array-operations":
  Load 1 ASDF system:
    array-operations
; Loading "array-operations"

(:ARRAY-OPERATIONS)

* (aops:generate #'identity 7 :position)
#(0 1 2 3 4 5 6)
~~~

Note that the nickname for `array-operations` is `aops`. The
`generate` function can also iterate over the array subscripts by
passing the key `:subscripts`. See the
[Array Operations manual on generate](https://lisp-stat.dev/docs/manuals/array-operations/#generate) for
more examples.

### Random numbers

To create an 3x3 array containing random numbers drawn from a uniform
distribution, `generate` can be used to call the CL
[random](http://clhs.lisp.se/Body/f_random.htm) function:

~~~lisp
* (aops:generate (lambda () (random 1.0)) '(3 3))
#2A((0.99292254 0.929777 0.93538976)
    (0.31522608 0.45167792 0.9411855)
    (0.96221936 0.9143338 0.21972346))
~~~

An array of Gaussian (normal) random numbers with mean of zero and
standard deviation of one, using the
[alexandria](https://common-lisp.net/project/alexandria/) package:

~~~lisp
* (ql:quickload :alexandria)
To load "alexandria":
  Load 1 ASDF system:
    alexandria
; Loading "alexandria"

(:ALEXANDRIA)

* (aops:generate #'alexandria:gaussian-random 4)
#(0.5522547885338768d0 -1.2564808468164517d0 0.9488161476129733d0
  -0.10372852118266523d0)
~~~

Note that this is not particularly efficient: It requires a function
call for each element, and although `gaussian-random` returns
two random numbers, only one of them is used.

For more efficient implementations, and a wider range of probability
distributions, there are packages available on Quicklisp. See
[CLiki](https://www.cliki.net/statistics) for a list.

## Accessing elements

To access the individual elements of an array there are the [aref](http://clhs.lisp.se/Body/f_aref.htm)
and [row-major-aref](http://clhs.lisp.se/Body/f_row_ma.htm#row-major-aref) functions.

The [aref](http://clhs.lisp.se/Body/f_aref.htm) function takes the same number of index arguments as the
array has dimensions. Indexing is from 0 and row-major as in C, but
not Fortran.

~~~lisp
* (defparameter *a* #(1 2 3 4))
*A*
* (aref *a* 0)
1
* (aref *a* 3)
4
* (defparameter *b* #2A((1 2 3) (4 5 6)))
*B*
* (aref *b* 1 0)
4
* (aref *b* 0 2)
3
~~~

The range of these indices can be found using
[array-dimensions](http://clhs.lisp.se/Body/f_ar_d_1.htm):

~~~
* (array-dimensions *a*)
(4)
* (array-dimensions *b*)
(2 3)
~~~

or the rank of the array can be found, and then the size of each
dimension queried:

~~~lisp
* (array-rank *a*)
1
* (array-dimension *a* 0)
4
* (array-rank *b*)
2
* (array-dimension *b* 0)
2
* (array-dimension *b* 1)
3
~~~

To loop over an array nested loops can be used, such as:

~~~lisp
* (defparameter a #2A((1 2 3) (4 5 6)))
A
* (destructuring-bind (n m) (array-dimensions a)
    (loop for i from 0 below n do
      (loop for j from 0 below m do
        (format t "a[~a ~a] = ~a~%" i j (aref a i j)))))

a[0 0] = 1
a[0 1] = 2
a[0 2] = 3
a[1 0] = 4
a[1 1] = 5
a[1 2] = 6
NIL
~~~

A utility macro which does this for multiple dimensions is `nested-loop`:

~~~lisp
(defmacro nested-loop (syms dimensions &body body)
  "Iterates over a multidimensional range of indices.

   SYMS must be a list of symbols, with the first symbol
   corresponding to the outermost loop.

   DIMENSIONS will be evaluated, and must be a list of
   dimension sizes, of the same length as SYMS.

   Example:
    (nested-loop (i j) '(10 20) (format t '~a ~a~%' i j))

  "
  (unless syms (return-from nested-loop `(progn ,@body))) ; No symbols

  ;; Generate gensyms for dimension sizes
  (let* ((rank (length syms))
         (syms-rev (reverse syms)) ; Reverse, since starting with innermost
         (dims-rev (loop for i from 0 below rank collecting (gensym))) ; innermost dimension first
         (result `(progn ,@body))) ; Start with innermost expression
    ;; Wrap previous result inside a loop for each dimension
    (loop for sym in syms-rev for dim in dims-rev do
         (unless (symbolp sym) (error "~S is not a symbol. First argument to nested-loop must be a list of symbols" sym))
         (setf result
               `(loop for ,sym from 0 below ,dim do
                     ,result)))
    ;; Add checking of rank and dimension types, and get dimensions into gensym list
    (let ((dims (gensym)))
      `(let ((,dims ,dimensions))
         (unless (= (length ,dims) ,rank) (error "Incorrect number of dimensions: Expected ~a but got ~a" ,rank (length ,dims)))
         (dolist (dim ,dims)
           (unless (integerp dim) (error "Dimensions must be integers: ~S" dim)))
         (destructuring-bind ,(reverse dims-rev) ,dims ; Dimensions reversed so that innermost is last
           ,result)))))
~~~

so that the contents of a 2D array can be printed using:

~~~lisp
* (defparameter a #2A((1 2 3) (4 5 6)))
A
* (nested-loop (i j) (array-dimensions a)
      (format t "a[~a ~a] = ~a~%" i j (aref a i j)))

a[0 0] = 1
a[0 1] = 2
a[0 2] = 3
a[1 0] = 4
a[1 1] = 5
a[1 2] = 6
NIL
~~~

[Note: This macro is available in [this fork](https://github.com/bendudson/array-operations) of array-operations, but
not Quicklisp]

### Row major indexing

In some cases, particularly element-wise operations, the number of
dimensions does not matter. To write code which is independent of the
number of dimensions, array element access can be done using a
single flattened index via
[row-major-aref](http://clhs.lisp.se/Body/f_row_ma.htm#row-major-aref).
The array size is given by [array-total-size](http://clhs.lisp.se/Body/f_ar_tot.htm), with the flattened
index starting at 0.

~~~lisp
* (defparameter a #2A((1 2 3) (4 5 6)))
A
* (array-total-size a)
6
* (loop for i from 0 below (array-total-size a) do
     (setf (row-major-aref a i) (+ 2.0 (row-major-aref a i))))
NIL
* a
#2A((3.0 4.0 5.0) (6.0 7.0 8.0))
~~~

### Infix syntax

The [cmu-infix](https://github.com/rigetticomputing/cmu-infix) library
provides some different syntax which can make mathematical expressions
easier to read:

~~~lisp
* (ql:quickload :cmu-infix)
To load "cmu-infix":
  Load 1 ASDF system:
    cmu-infix
; Loading "cmu-infix"

(:CMU-INFIX)

* (named-readtables:in-readtable cmu-infix:syntax)
(("COMMON-LISP-USER" . #<NAMED-READTABLE CMU-INFIX:SYNTAX {10030158B3}>)
 ...)

* (defparameter arr (make-array '(3 2) :initial-element 1.0))
ARR

* #i(arr[0 1] = 2.0)
2.0

* arr
#2A((1.0 2.0) (1.0 1.0) (1.0 1.0))
~~~

A matrix-matrix multiply operation can be implemented as:

~~~lisp
(let ((A #2A((1 2) (3 4)))
      (B #2A((5 6) (7 8)))
      (result (make-array '(2 2) :initial-element 0.0)))

     (loop for i from 0 to 1 do
           (loop for j from 0 to 1 do
                 (loop for k from 0 to 1 do
                       #i(result[i j] += A[i k] * B[k j]))))
      result)
~~~

See the section below on linear algebra, for alternative
matrix-multiply implementations.

## Element-wise operations

To multiply two arrays of numbers of the same size, pass a function
to `each` in the [array-operations](https://github.com/Lisp-Stat/array-operations) library:

~~~lisp
* (aops:each #'* #(1 2 3) #(2 3 4))
#(2 6 12)
~~~

For improved efficiency there is the `aops:each*` function, which
takes a type as first argument to specialise the result array.

To add a constant to all elements of an array:

~~~lisp
* (defparameter *a* #(1 2 3 4))
*A*
* (aops:each (lambda (it) (+ 42 it)) *a*)
#(43 44 45 46)
* *a*
#(1 2 3 4)
~~~

Note that `each` is not destructive, but makes a new array.
All arguments to `each` must be arrays of the same size,
so `(aops:each #'+ 42 *a*)` is not valid.

### Vectorising expressions

An alternative approach to the `each` function above, is to use a
macro to iterate over all elements of an array:

~~~lisp
(defmacro vectorize (variables &body body)
  ;; Check that variables is a list of only symbols
  (dolist (var variables)
    (if (not (symbolp var))
        (error "~S is not a symbol" var)))

    ;; Get the size of the first variable, and create a new array
    ;; of the same type for the result
    `(let ((size (array-total-size ,(first variables)))  ; Total array size (same for all variables)
           (result (make-array (array-dimensions ,(first variables)) ; Returned array
                               :element-type (array-element-type ,(first variables)))))
       ;; Check that all variables have the same sizeo
       ,@(mapcar (lambda (var) `(if (not (equal (array-dimensions ,(first variables))
                                                (array-dimensions ,var)))
                                    (error "~S and ~S have different dimensions" ',(first variables) ',var)))
              (rest variables))

       (dotimes (indx size)
         ;; Locally redefine variables to be scalars at a given index
         (let ,(mapcar (lambda (var) (list var `(row-major-aref ,var indx))) variables)
           ;; User-supplied function body now evaluated for each index in turn
           (setf (row-major-aref result indx) (progn ,@body))))
       result))
~~~

[Note: Expanded versions of this macro are available in [this
fork](https://github.com/bendudson/array-operations) of array-operations, but
not Quicklisp]

This can be used as:

~~~lisp
* (defparameter *a* #(1 2 3 4))
*A*
* (vectorize (*a*) (* 2 *a*))
#(2 4 6 8)
~~~

Inside the body of the expression (second form in `vectorize`
expression) the symbol `*a*` is bound to a single element.
This means that the built-in mathematical functions can be used:

~~~lisp
* (defparameter a #(1 2 3 4))
A
* (defparameter b #(2 3 4 5))
B
* (vectorize (a b) (* a (sin b)))
#(0.9092974 0.28224 -2.2704074 -3.8356972)
~~~

and combined with `cmu-infix`:

~~~lisp
* (vectorize (a b) #i(a * sin(b)) )
#(0.9092974 0.28224 -2.2704074 -3.8356972)
~~~

### Calling BLAS

Several packages provide wrappers around BLAS, for fast matrix manipulation.

The [lla](https://github.com/tpapp/lla) package in quicklisp includes
calls to some functions:

#### Scale an array

scaling by a constant factor:

~~~lisp
* (defparameter a #(1 2 3))
* (lla:scal! 2.0 a)
* a
#(2.0d0 4.0d0 6.0d0)
~~~

#### AXPY

This calculates `a * x + y` where `a` is a constant, `x` and `y` are arrays.
The `lla:axpy!` function is destructive, modifying the last argument (`y`).

~~~lisp
* (defparameter x #(1 2 3))
A
* (defparameter y #(2 3 4))
B
* (lla:axpy! 0.5 x y)
#(2.5d0 4.0d0 5.5d0)
* x
#(1.0d0 2.0d0 3.0d0)
* y
#(2.5d0 4.0d0 5.5d0)
~~~

If the `y` array is complex, then this operation calls the complex
number versions of these operators:

~~~lisp
* (defparameter x #(1 2 3))
* (defparameter y (make-array 3 :element-type '(complex double-float)
                                :initial-element #C(1d0 1d0)))
* y
#(#C(1.0d0 1.0d0) #C(1.0d0 1.0d0) #C(1.0d0 1.0d0))

* (lla:axpy! #C(0.5 0.5) a b)
#(#C(1.5d0 1.5d0) #C(2.0d0 2.0d0) #C(2.5d0 2.5d0))
~~~

#### Dot product

The dot product of two vectors:

~~~lisp
* (defparameter x #(1 2 3))
* (defparameter y #(2 3 4))
* (lla:dot x y)
20.0d0
~~~

### Reductions

The [reduce](http://clhs.lisp.se/Body/f_reduce.htm) function operates
on sequences, including vectors (1D arrays), but not on
multidimensional arrays.
To get around this, multidimensional arrays can be displaced to create
a 1D vector.
Displaced arrays share storage with the original array, so this is a
fast operation which does not require copying data:

~~~lisp
* (defparameter a #2A((1 2) (3 4)))
A
* (reduce #'max (make-array (array-total-size a) :displaced-to a))
4
~~~

The `array-operations` package contains `flatten`, which returns a
displaced array i.e doesn't copy data:

~~~lisp
* (reduce #'max (aops:flatten a))
~~~

An SBCL extension,
[array-storage-vector](http://www.sbcl.org/manual/#index-array_002dstorage_002dvector)
provides an efficient but not portable way to achieve the same thing:

~~~lisp
* (reduce #'max (array-storage-vector a))
4
~~~

More complex reductions are sometimes needed, for example finding the
maximum absolute difference between two arrays. Using the above
methods we could do:

~~~lisp
* (defparameter a #2A((1 2) (3 4)))
A
* (defparameter b #2A((1 3) (5 4)))
B
* (reduce #'max (aops:flatten
                  (aops:each
                    (lambda (a b) (abs (- a b))) a b)))
2
~~~

This involves allocating an array to hold the intermediate result,
which for large arrays could be inefficient. Similarly to `vectorize`
defined above, a macro which does not allocate can be defined as:

~~~lisp
(defmacro vectorize-reduce (fn variables &body body)
  "Performs a reduction using FN over all elements in a vectorized expression
   on array VARIABLES.

   VARIABLES must be a list of symbols bound to arrays.
   Each array must have the same dimensions. These are
   checked at compile and run-time respectively.
  "
  ;; Check that variables is a list of only symbols
  (dolist (var variables)
    (if (not (symbolp var))
        (error "~S is not a symbol" var)))

  (let ((size (gensym)) ; Total array size (same for all variables)
        (result (gensym)) ; Returned value
        (indx (gensym)))  ; Index inside loop from 0 to size

    ;; Get the size of the first variable
    `(let ((,size (array-total-size ,(first variables))))
       ;; Check that all variables have the same size
       ,@(mapcar (lambda (var) `(if (not (equal (array-dimensions ,(first variables))
                                                (array-dimensions ,var)))
                                    (error "~S and ~S have different dimensions" ',(first variables) ',var)))
              (rest variables))

       ;; Apply FN with the first two elements (or fewer if size < 2)
       (let ((,result (apply ,fn (loop for ,indx below (min ,size 2) collecting
                                      (let ,(map 'list (lambda (var) (list var `(row-major-aref ,var ,indx))) variables)
                                        (progn ,@body))))))

         ;; Loop over the remaining indices
         (loop for ,indx from 2 below ,size do
            ;; Locally redefine variables to be scalars at a given index
              (let ,(mapcar (lambda (var) (list var `(row-major-aref ,var ,indx))) variables)
                ;; User-supplied function body now evaluated for each index in turn
                (setf ,result (funcall ,fn ,result (progn ,@body)))))
         ,result))))

~~~

[Note: This macro is available in [this fork](https://github.com/bendudson/array-operations)
of array-operations, but not Quicklisp]

Using this macro, the maximum value in an array A (of any shape) is:

~~~lisp
* (vectorize-reduce #'max (a) a)
~~~

The maximum absolute difference between two arrays A and B, of any
shape as long as they have the same shape, is:

~~~lisp
* (vectorize-reduce #'max (a b) (abs (- a b)))
~~~

## Linear algebra

Several packages provide bindings to BLAS and LAPACK libraries,
including:

- [lla](https://github.com/tpapp/lla)
- [MAGICL](https://github.com/rigetticomputing/magicl)

A longer list of available packages is on [CLiki's linear algebra page](https://www.cliki.net/linear%20algebra).

In the examples below the lla package is loaded:

~~~lisp
* (ql:quickload :lla)

To load "lla":
  Load 1 ASDF system:
    lla
; Loading "lla"
.
(:LLA)
~~~

### Matrix multiplication

The [lla](https://github.com/tpapp/lla) function `mm` performs
vector-vector, matrix-vector and matrix-matrix
multiplication.

#### Vector dot product

Note that one vector is treated as a row vector, and the other as
column:

~~~lisp
* (lla:mm #(1 2 3) #(2 3 4))
20
~~~

#### Matrix-vector product

~~~lisp
* (lla:mm #2A((1 1 1) (2 2 2) (3 3 3))  #(2 3 4))
#(9.0d0 18.0d0 27.0d0)
~~~

which has performed the sum over `j` of `A[i j] * x[j]`

#### Matrix-matrix multiply

~~~lisp
* (lla:mm #2A((1 2 3) (1 2 3) (1 2 3))  #2A((2 3 4) (2 3 4) (2 3 4)))
#2A((12.0d0 18.0d0 24.0d0) (12.0d0 18.0d0 24.0d0) (12.0d0 18.0d0 24.0d0))
~~~

which summed over `j` in `A[i j] * B[j k]`

Note that the type of the returned arrays are simple arrays,
specialised to element type `double-float`

~~~lisp
* (type-of (lla:mm #2A((1 0 0) (0 1 0) (0 0 1)) #(1 2 3)))
(SIMPLE-ARRAY DOUBLE-FLOAT (3))
~~~

#### Outer product

The [array-operations](https://github.com/Lisp-Stat/array-operations)
package contains a generalised [outer product](https://en.wikipedia.org/wiki/Outer_product)
function:

~~~lisp
* (ql:quickload :array-operations)
To load "array-operations":
  Load 1 ASDF system:
    array-operations
; Loading "array-operations"

(:ARRAY-OPERATIONS)
* (aops:outer #'* #(1 2 3) #(2 3 4))
#2A((2 3 4) (4 6 8) (6 9 12))
~~~

which has created a new 2D array `A[i j] = B[i] * C[j]`. This `outer`
function can take an arbitrary number of inputs, and inputs with
multiple dimensions.

### Matrix inverse

The direct inverse of a dense matrix can be calculated with `invert`

~~~lisp
* (lla:invert #2A((1 0 0) (0 1 0) (0 0 1)))
#2A((1.0d0 0.0d0 -0.0d0) (0.0d0 1.0d0 -0.0d0) (0.0d0 0.0d0 1.0d0))
~~~

e.g

~~~lisp
* (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
A
* (defparameter b (lla:invert a))
B
* (lla:mm a b)
#2A((1.0d0 2.220446049250313d-16 0.0d0)
    (0.0d0 1.0d0 0.0d0)
    (0.0d0 1.1102230246251565d-16 0.9999999999999998d0))
~~~

Calculating the direct inverse is generally not advisable,
particularly for large matrices. Instead the
[LU decomposition](https://en.wikipedia.org/wiki/LU_decomposition) can be
calculated and used for multiple inversions.

~~~lisp
* (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
A
* (defparameter b (lla:mm a #(1 2 3)))
B
* (lla:solve (lla:lu a) b)
#(1.0d0 2.0d0 3.0d0)
~~~

### Singular value decomposition

The `svd` function calculates the [singular value decomposition](https://en.wikipedia.org/wiki/Singular-value_decomposition)
of a given matrix, returning an object with slots for the three returned
matrices:

~~~lisp
* (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
A
* (defparameter a-svd (lla:svd a))
A-SVD
* a-svd
#S(LLA:SVD
   :U #2A((-0.6494608633564334d0 0.7205486773948702d0 0.24292013188045855d0)
          (-0.3744175632000917d0 -0.5810891192666799d0 0.7225973455785591d0)
          (-0.6618248071322363d0 -0.3783451320875919d0 -0.6471807210432038d0))
   :D #S(CL-NUM-UTILS.MATRIX:DIAGONAL-MATRIX
         :ELEMENTS #(5.593122609997059d0 1.2364443401235103d0
                     0.43380279311714376d0))
   :VT #2A((-0.2344460799312531d0 -0.7211054639318696d0 -0.6519524104506949d0)
           (0.2767642134809678d0 -0.6924017945853318d0 0.6663192365460215d0)
           (-0.9318994611765425d0 -0.02422116311440764d0 0.3619070730398283d0)))
~~~

The diagonal matrix (singular values) and vectors can be accessed with
functions:

~~~lisp
(lla:svd-u a-svd)
#2A((-0.6494608633564334d0 0.7205486773948702d0 0.24292013188045855d0)
    (-0.3744175632000917d0 -0.5810891192666799d0 0.7225973455785591d0)
    (-0.6618248071322363d0 -0.3783451320875919d0 -0.6471807210432038d0))

* (lla:svd-d a-svd)
#S(CL-NUM-UTILS.MATRIX:DIAGONAL-MATRIX
   :ELEMENTS #(5.593122609997059d0 1.2364443401235103d0 0.43380279311714376d0))

* (lla:svd-vt a-svd)
#2A((-0.2344460799312531d0 -0.7211054639318696d0 -0.6519524104506949d0)
    (0.2767642134809678d0 -0.6924017945853318d0 0.6663192365460215d0)
    (-0.9318994611765425d0 -0.02422116311440764d0 0.3619070730398283d0))
~~~


## Matlisp

The [Matlisp](https://github.com/bharath1097/matlisp/) scientific
computation library provides high performance operations on arrays,
including wrappers around BLAS and LAPACK functions.
It can be loaded using quicklisp:

~~~lisp
* (ql:quickload :matlisp)
~~~

The nickname for `matlisp` is `m`. To avoid typing `matlisp:` or
`m:` in front of each symbol, you can define your own package which
uses matlisp
(See the [PCL section on packages](http://www.gigamonkeys.com/book/programming-in-the-large-packages-and-symbols.html)):

~~~lisp
* (defpackage :my-new-code
     (:use :common-lisp :matlisp))
#<PACKAGE "MY-NEW-CODE">

* (in-package :my-new-code)
~~~

and to use the `#i` infix reader (note the same name as for
`cmu-infix`), run:

~~~lisp
* (named-readtables:in-readtable :infix-dispatch-table)
~~~

### Creating tensors

~~~lisp
* (matlisp:zeros '(2 2))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
  0.000    0.000
  0.000    0.000
>
~~~

Note that by default matrix storage types are `double-float`.
To create a complex array using `zeros`, `ones` and `eye`, specify the type:

~~~lisp
* (matlisp:zeros '(2 2) '((complex double-float)))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: (COMPLEX DOUBLE-FLOAT)>| #(2 2)
  0.000    0.000
  0.000    0.000
>
~~~

As well as `zeros` and `ones` there is `eye` which creates an identity
matrix:

~~~lisp
* (matlisp:eye '(3 3) '((complex double-float)))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: (COMPLEX DOUBLE-FLOAT)>| #(3 3)
  1.000    0.000    0.000
  0.000    1.000    0.000
  0.000    0.000    1.000
>
~~~

#### Ranges

To generate 1D arrays there are the `range` and `linspace` functions:

~~~lisp
* (matlisp:range 1 10)
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(9)
 1   2   3   4   5   6   7   8   9
>
~~~

The `range` function rounds down it's final argument to an integer:

~~~lisp
* (matlisp:range 1 -3.5)
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: SINGLE-FLOAT>| #(5)
 1.000   0.000   -1.000  -2.000  -3.000
>
* (matlisp:range 1 3.3)
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: SINGLE-FLOAT>| #(3)
 1.000   2.000   3.000
>
~~~

Linspace is a bit more general, and the values returned include the
end point.

~~~lisp
* (matlisp:linspace 1 10)
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(10)
 1   2   3   4   5   6   7   8   9   10
>
~~~

~~~lisp
* (matlisp:linspace 0 (* 2 pi) 5)
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(5)
 0.000   1.571   3.142   4.712   6.283
>
~~~

Currently `linspace` requires real inputs, and doesn't work with complex numbers.

#### Random numbers

~~~lisp
* (matlisp:random-uniform '(2 2))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
  0.7287       0.9480
  2.6703E-2    0.1834
>
~~~

~~~lisp
(matlisp:random-normal '(2 2))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
  0.3536    -1.291
 -0.3877    -1.371
>
~~~

There are functions for other distributions, including
`random-exponential`, `random-beta`, `random-gamma` and
`random-pareto`.

#### Reader macros

The `#d` and `#e` reader macros provide a way to create `double-float`
and `single-float` tensors:

~~~lisp
* #d[1,2,3]
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(3)
 1.000   2.000   3.000
>

* #d[[1,2,3],[4,5,6]]
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
  1.000    2.000    3.000
  4.000    5.000    6.000
>
~~~

Note that the comma separators are needed.

#### Tensors from arrays

Common lisp arrays can be converted to Matlisp tensors by copying:

~~~lisp
* (copy #2A((1 2 3)
            (4 5 6))
        '#.(tensor 'double-float))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
  1.000    2.000    3.000
  4.000    5.000    6.000
>
~~~

Instances of the `tensor` class can also be created, specifying the
dimensions. The internal storage of `tensor` objects is a 1D array
(`simple-vector`) in a slot `store`.

For example, to create a `double-float` type tensor:

~~~lisp
(make-instance (tensor 'double-float)
    :dimensions  (coerce '(2) '(simple-array index-type (*)))
    :store (make-array 2 :element-type 'double-float))
~~~

#### Arrays from tensors

The array store can be accessed using slots:

~~~lisp
* (defparameter vec (m:range 0 5))
* vec
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(5)
 0   1   2   3   4
>
* (slot-value vec 'm:store)
#(0 1 2 3 4)
~~~

Multidimensional tensors are also stored in 1D arrays, and are stored
in column-major order rather than the row-major ordering used for
common lisp arrays. A displaced array will therefore be
transposed.

The contents of a tensor can be copied into an array

~~~lisp
* (let ((tens (m:ones '(2 3))))
    (m:copy tens 'array))
#2A((1.0d0 1.0d0 1.0d0) (1.0d0 1.0d0 1.0d0))
~~~

or a list:

~~~lisp
* (m:copy (m:ones '(2 3)) 'cons)
((1.0d0 1.0d0 1.0d0) (1.0d0 1.0d0 1.0d0))
~~~

### Element access

The `ref` function is the equivalent of `aref` for standard CL
arrays, and is also setf-able:

~~~lisp
* (defparameter a (matlisp:ones '(2 3)))

* (setf (ref a 1 1) 2.0)
2.0d0
* a
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
  1.000    1.000    1.000
  1.000    2.000    1.000
>
~~~

### Element-wise operations

The `matlisp-user` package, loaded when `matlisp` is loaded, contains
functions for operating element-wise on tensors.

~~~lisp
* (matlisp-user:* 2 (ones '(2 3)))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
  2.000    2.000    2.000
  2.000    2.000    2.000
>
~~~

This includes arithmetic operators '+', '-', '*', '/' and 'expt', but
also `sqrt`,`sin`,`cos`,`tan`, hyperbolic functions, and their inverses.
The `#i` reader macro recognises many of these, and uses the
`matlisp-user` functions:

~~~lisp
* (let ((a (ones '(2 2)))
        (b (random-normal '(2 2))))
     #i( 2 * a + b ))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
  0.9684    3.250
  1.593     1.508
>

* (let ((a (ones '(2 2)))
        (b (random-normal '(2 2))))
     (macroexpand-1 '#i( 2 * a + b )))
(MATLISP-USER:+ (MATLISP-USER:* 2 A) B)
~~~
